{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Варахобова Анна Андреевна, @h.varakhobava."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Обзор библиотеки статистического моделирования Statsmodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В этом тьюториале будет рассмотрена библиотека statsmodels и некоторые популярные ее применения.\n",
    "\n",
    "<a href=\"https://www.statsmodels.org\">statsmodels</a> - это библиотека для статистического моделирования, содержащая всевозможные статистические модели, множество статистических тестов, средств для информативных графиков и др. Распространяется под лицензией New BSD (modified BSD).\n",
    "\n",
    "Этот пакет покажется знакомым тем, кто работал со статистическим моделированием на R, т.к. позволяет пользоватся строковыми формулами для описания линейной модели, и содержит похожие методы и классы, а также позволяет пользоваться наборами данных из R.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Линейная регрессия"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, как решить задачу линейной регрессии с помощью statsmodels. В пакете statsmodels.api.datasets представлено некоторое количество датасетов, включая датасеты, доступные в R. Воспользуемся одним из них. \n",
    "\n",
    "Во встроенных данных не из R присутствуют свойства endog и exog - целевой и все остальные признаки соответственно, а сам массив лежит в свойстве data. В даных из R endog и exog свойств нет, а в data лежит pandas DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.api as sms\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "from sklearn.linear_model import LinearRegression, LogisticRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.get_rdataset(\"Guerry\", \"HistData\").data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['Lottery'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для обычной линейной регрессии используем OLS - ordinary least squares модель."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сначала подготовим данные: удалим пропуски и закодируем категориальный признак Region, а затем обучим модель и выведем саммари."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = data[['Lottery', 'Literacy', 'Wealth', 'Region']].dropna()\n",
    "\n",
    "df_dumm = pd.get_dummies(df, columns=['Region'], drop_first=True)\n",
    "\n",
    "# обучим, добавив к предикторам смещение\n",
    "ols = sm.OLS(df_dumm['Lottery'], sm.add_constant(df_dumm.drop(['Lottery'], axis=1)))\n",
    "res = ols.fit()\n",
    "\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В саммари можно сразу без лишних телодвижений посмотреть статистические характеристики модели: коэффициенты модели, значимы ли они, их доверительные интервалы. Также выводится множество других параметров: метрики качества, значения статистических тестов.\n",
    "\n",
    "Видим, что модель имеет R2 метрику, равную 0.3, так себе модель :) Высокие p-value коэффициентов говорят о том же."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Формулы\n",
    "Проделаем то же самое, но с использованием формул, а также посмотрим, какие \"чудеса\" можно вытворять с их помощью.\n",
    "\n",
    "Говорим модели, что хотим предсказать Lottery как линейную комбинацию остальных перечисленных параметров\n",
    "добавляем вконце -1, чтобы убрать смещение, т.е. чтобы y = a1x1 + .. + anxn\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_f = smf.ols(formula='Lottery ~ Literacy + Wealth + Region', data=df).fit()\n",
    "res_f.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, тот же результат получен более удобно, без ручных преобразований признаков. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С помощью формул также можно катеригоризировать признаки, удалять сонстанту (смещение), конструировать произведение признаков, а также применять к ним функции. Для экономии пространства будем выводить только параметры модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# категоризация\n",
    "res_f = smf.ols(formula='Lottery ~ Literacy + Wealth + C(Region)', data=df).fit()\n",
    "res_f.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# удаляем смещение\n",
    "res_f = smf.ols(formula='Lottery ~ Literacy + Wealth + C(Region) - 1', data=df).fit()\n",
    "res_f.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# заменим два признака на их произведение\n",
    "res_f = smf.ols(formula='Lottery ~ Literacy:Wealth + C(Region) - 1', data=df).fit()\n",
    "res_f.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# добавим произведение двух признаков, но оставим исходные\n",
    "res_f = smf.ols(formula='Lottery ~ Literacy*Wealth + C(Region) - 1', data=df).fit()\n",
    "res_f.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# применим функцию к признаку\n",
    "res_f = smf.ols(formula='Lottery ~ np.log(Literacy) + Wealth', data=df).fit()\n",
    "res_f.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Но самая полезная киллер-фича формул - возможность использовать их в моделях, не поддреживающих такую запись. Для этого нам понадобится библиотека patsy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import patsy\n",
    "\n",
    "f = 'Lottery ~ Literacy * Wealth + C(Region) - 1'\n",
    "\n",
    "y, X = patsy.dmatrices(f, df, return_type='dataframe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Робастная регрессия\n",
    "\n",
    "Иногда в данных присутсвуют выбросы, которые влияют на качество модели. В таких случаях можно использовать линейную модель, устойчивую к ним, - RLM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "\n",
    "new_df = pd.DataFrame()\n",
    "new_df['x1'] = np.linspace(0, 20, nsample)\n",
    "new_df['x2'] = (new_df['x1'] - 5)**2\n",
    "\n",
    "sig = 0.3   # маленькая вариация, больше разницы между моделями OLS и RLM\n",
    "beta = [0.5, -0.0]\n",
    "kx = np.dot(new_df.values, beta)\n",
    "y= kx + sig*1. * np.random.normal(size=nsample)\n",
    "y_true = kx.copy() # значения без выбросов\n",
    "y[[38,40,42,45,49]] -= 5  # выбросы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df['y'] = y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols_res = smf.ols('y ~ x1 + x2', data=new_df).fit()\n",
    "ols_res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обучим RLM модель"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_rlm = smf.rlm('y ~ x1 + x2', data=new_df).fit()\n",
    "res_rlm.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df['fitted_ols'] = ols_res.fittedvalues\n",
    "new_df['fitted_rlm'] = res_rlm.fittedvalues\n",
    "\n",
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot(new_df.x1.values, y, 'o', c='b', label='y')\n",
    "plt.plot(new_df.x1.values, y_true, c='black', label='True')\n",
    "plt.plot(new_df.x1, new_df.fitted_ols, c='r', label='OLS')\n",
    "plt.plot(new_df.x1, new_df.fitted_rlm, c='g', label='RLM')\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, RLM справляется с выбросами гораздо лучше OLS. В качестве робастной нормы можно передавать LeastSquares, HuberT, RamsayE, AndrewWave, TrimmedMean, Hampel, and TukeyBiweight из пакета sm.robust.norms (подробнее <a href=\"https://www.statsmodels.org/stable/rlm.html\">тут</a> )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GLM и компания"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Линейная регрессия накладывает на данные некоторые ограничения, в частности, нормальность распределения ошибок и целевой переменной, а также линейность предсказаний. Но что если это не так? Можно долго мучаться с преобразованием данных, а можно посмотреть в сторону GLM - <a href='https://en.wikipedia.org/wiki/Generalized_linear_model'>generalized linear model</a>. GLM часто применяется в страховании для моделирования частоты заявлений об ущербе и среднего размера выплат, чтобы рассчитать размер взноса за полис, а также в других случаях, когда целевая переменная распределена не по нормальному закону.\n",
    "<br>\n",
    "<br>\n",
    "Рассмотрим применение модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = sm.datasets.star98.load()\n",
    "print(sm.datasets.star98.NOTE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.data[:1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.endog[:3], ds.exog[:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обучим GLM, указав распределениe sm.families.Binomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = sm.GLM(ds.endog, ds.exog, family=sm.families.Binomial()).fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на остатки модели, для этого построим гистограмму распределения дисперсии остатков и построим Q-Q plot с помощью statsmodels.graphics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "plt.hist(res.resid_deviance.copy(), bins=25)\n",
    "plt.title('Histogram of standardized deviance residuals');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels import graphics\n",
    "\n",
    "graphics.gofplots.qqplot(resid, line='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что остатки распределены нормально. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В statsmodels также реализованы расширения GLM: **GEE**, **MixedLM** - для кластеризованных данных, не скореллированных между кластерами, но внутри кластеров, и с наличием случайных корреляций неизвестной природы. Подробнее <a href=\"https://www.statsmodels.org/stable/gee.html\">тут</a> и <a href = \"https://www.statsmodels.org/stable/mixed_linear.html\">тут</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Логистическая регрессия\n",
    "\n",
    "Посмотрим, как применить statsmodels к задаче классификации:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spector_data = sm.datasets.spector.load()\n",
    "print(sm.datasets.spector.NOTE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spector_data.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### если при запуске выдается ошибка \"module 'scipy.stats' has no attribute 'chisqprob'\", \n",
    "### раскомментируйте строки ниже\n",
    "# from scipy import stats\n",
    "# stats.chisqprob = lambda chisq, df: stats.chi2.sf(chisq, df)\n",
    "\n",
    "res_logit = smf.logit('GRADE ~ TUCE + PSI + GPA', data=spector_data.data).fit()\n",
    "res_logit.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(spector_data.data['PSI'], spector_data.data['GRADE'], 'o', c='b')\n",
    "plt.plot(spector_data.data['PSI'], res_logit.fittedvalues, 'o', c='g')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EDA\n",
    "В модуле statsmodels.grphics есть много других возможностей для визуального анализа данных, посмотрим на некоторые из них."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.exog.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.graphics import boxplots, correlation, regressionplots\n",
    "\n",
    "boxplots.violinplot(ds.exog[:, :5], positions=[0, 1, 2, 3, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = np.corrcoef(ds.exog[:5])\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "correlation.plot_corr(cm, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Мощный тул для отображения результатов регрессии:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure(figsize=(10, 10))\n",
    "regressionplots.plot_regress_exog(ols_res, 1, fig=f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Полезные функции"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В библиотеке много полезных функций, сосредоточены они в модуле statsmodels.tools.tools."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Категоризация:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tools import tools\n",
    "\n",
    "tools.categorical(df.Region.values)[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Добавление смещения к признакам, полезно, если исплользуем модель без формул:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tools.add_constant(df).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Метрики точности предсказания:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tools import eval_measures\n",
    "\n",
    "ols_rmse = eval_measures.rmse(y, ols_res.fittedvalues)\n",
    "ols_mse = eval_measures.mse(y, ols_res.fittedvalues)\n",
    "print('OLS metrics: {}, {}'.format(ols_rmse, ols_mse))\n",
    "\n",
    "rlm_rmse = eval_measures.rmse(y_true, res_rlm.fittedvalues)\n",
    "rlm_mse = eval_measures.mse(y_true, res_rlm.fittedvalues)\n",
    "print('RLM metrics: {}, {}'.format(rlm_rmse, rlm_mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Заключение\n",
    "\n",
    "В этом тьюториале мы рассмотрели самые расспространенные приемы моделирования и другие полезные инструменты statsmodels. В библиотеке еще много всего: ANOVA, непараметрические методы, методы прогнозирования временных рядов, статистические тесты.\n",
    "\n",
    "В рассмотренных случаях библиотека имеет ряд преимуществ перед slkearn:\n",
    "* лаконичность\n",
    "* развернутое саммари по модели из коробки\n",
    "* формулы\n",
    "* множество статистических интсрументов и проверок\n",
    "* минималистичность в целом\n",
    "\n",
    "Так что если вам не нужно что-то сложное типа градиентного бустинга, statsmodels отлично подойдет\n",
    "\n",
    "Спасибо за внимание!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
